Directed-loop Monte Carlo simulations of vertex models 
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We show how the directed-loop Monte Carlo algorithm can be applied to study vertex models. 
The algorithm is employed to calculate the arrow polarization in the six-vertex model with the 
domain wall boundary conditions (DWBC). The model exhibits spatially separated ordered and 
"disordered" regions. We show how the boundary between these regions depends on parameters 
of the model. We give some predictions on the behavior of the polarization in the thermodynamic 
limit and discuss the relation to the Arctic Circle theorem. 
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I. INTRODUCTION 



Vertex models have a long and distinguished history 
in physics. Their fame is intimately connected to the 
concept of integrability, and the exact solutions of the 
six-vertex pj and the eight-vertex models with peri- 
odic boundary conditions (PBC) are indeed milestones in 
physics. Despite being exactly solvable, there are ques- 
tions about these models that cannot easily be answered. 
An example is the influence of boundary conditions on 
correlation functions. While boundary conditions are not 
normally important in the thermodynamic limit, they 
have a profound influence on the vertex models. Exact 
studies, made for the six-vertex model with the domain 
wall boundary conditions (DWBC) yj show this in par- 
ticular. These studies were restricted to certain points 
in the phase diagram, and involve rather sophisticated 
mathematical methods. It is thus appropriate to com- 
plement them with Monte Carlo simulations. 

The purpose of this article is to demonstrate that the 
directed-loop Monte Carlo algorithm developed for quan- 
tum spin systems Q can be used as an effective tool to 
study vertex models. The discussion of the algorithm will 
be kept general, but when demonstrating its use we will 
focus on the six-vertex model with the DWBC, a model 
which is difficult to simulate using other known Monte 
Carlo algorithms. 



II. MONTE CARLO ALGORITHM 

In a vertex model, each vertex have edges with an 
Ising-like variable, an arrow, that points either away from 
or into the vertex. The arrangement of arrows around 
the vertex determines the vertex weight. Two vertices 
are joined by their common edge, sharing the arrow on 
the edge. In general there are no restrictions on which 
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FIG. 1: Illustration of the directed-loop algorithm. Vertex 
edges are drawn with two arrows allowing the discontinuity 
at the head and tail of the loop to be shown. The thick line 
shows the loop path along which the arrows has been nipped. 
The loop closes when the loop head (thick arrow) hits the 
loop tail (vertical bar). 



vertices are joined, however for traditional vertex models 
nearest-neighbor vertices are joined together. The Monte 
Carlo algorithm discussed here always flips two (or zero) 
arrows on a vertex, thus it is limited to models where 
an even number of arrows are pointing away from each 
vertex. Most vertex models of interest obey this rule. 

In visualizing the directed-loop Monte Carlo algorithm, 
originally developed for quantum systems 0] , it is helpful 
to cut every edge into two pieces, each piece having an 
arrow belonging to a specific vertex, Fig. ^ For a valid 
vertex configuration the arrows on the two parts of an 
edge must have the same orientation. The directed-loop 
algorithm is as follows: Pick a random vertex v\ and a 
random edge belonging to that vertex. Based on these 
choices select in a probabilistic manner another edge be- 
longing to v\ and name that the out-edge. Then flip the 
arrows on both the part of the in-edge and the part of the 
out-edge belonging to v\. This introduces two disconti- 
nuities in the arrow configurations on the edges, one on 
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FIG. 2: The vertices of the eight- vertex model and their sta- 
tistical weights. 



the starting in-edge and another one on the out-edge. 
The new configuration is thus not an allowed vertex con- 
figuration. To repair this, the out-edge discontinuity is 
moved by repeating the procedure on the vertex con- 
nected to the out-edge v 2 , this time using the out-edge of 
vi and the in-edge on v 2 . The process is stopped when 
the out-edge selected is the starting edge, thus healing all 
discontinuities. In this way arrows are flipped as a loop 
is constructed, and a new allowed vertex configuration is 
arrived at when the loop closes. 

In order to determine the probabilities for selecting 
out-edges and to see how detailed balance is satisfied one 
needs to consider also the probability for the reverse up- 
date. The reverse update consists of traversing the same 
loop in the opposite direction while flipping arrows back. 
As is explained in detail in Ref. P|, detailed balance is 
satisfied for the whole loop construction, if detailed bal- 
ance is satisfied in each edge selecting step, for which the 
criterion is as follows: Let w be the weight of the ver- 
tex v before edge-flips, then the probability P(v, i —> o) 
for exiting at the out-edge o, given that the in-edge is i, 
should satisfy 
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where w' is the weight of the vertex v obtained by flip- 
ping the arrows on edges i and o belonging to the vertex 
v. Notice that P(v',o — > i), on the right hand side, de- 
scribes an edge-selecting step in the reverse update pro- 
cess where the loop is traversed in the opposite direction 
to that described on the left hand side. The loop con- 
struction should not terminate in the edge-selecting step, 
thus 
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where the sum is taken over all possible out-edges, in- 
cluding the in-edge i. 

This algorithm resembles closely the ice model algo- 
rithm invented by Rahman and Stillinger [{| , generalized 
to arbitrary couplings by Barkema and Newman [(|. In 
fact, at the point in parameter space where all vertex 
weights are equal our algorithm is identical to the long- 
loop version of the ice model algorithm. However away 



from this point, Barkema and Newman's algorithm in- 
volves accepting or rejecting the loop after it has been 
constructed. The directed-loop algorithm has no such 
accept/reject step. A comparison of integrated autocor- 
relation times for the directed-loop algorithm and the 
short-loop algorithm of Barkema and Newman are shown 
in Fig. |21 The autocorrelation times are measured in 
units of lattice sweeps. One lattice sweep corresponds to 
a number of elementary loop moves such that on average 
each vertex on the lattice have been visited twice. In 
defining visited we include parts of the loop where the 
loop bounces off a vertex (relevant for the directed-loop 
algorithm) and the neck part of short-loops. Neither the 
bounces nor the short-loop-necks contribute to changes 
in the vertex configuration. However they are intrinsic 
parts of the algorithms and requires computer effort, and 
should therefore be accounted for. 

The upper panel of Fig. [31 shows integrated autocor- 
relation times of the observable counting the number of 
c-type vertices in each configuration. This observables 
was chosen to compare with the performance results in 
Ref. While the integrated autocorrelation times are 
larger for the short-loop algorithm the scaling with sys- 
tem size appears to be equal for both algorithms. The 
lower panel shows integrated autocorrelation times for 
the total arrow-polarization in the y-direction. These 
scales much worse for the short-loop algorithm than for 
the directed-loop algorithm. This is to be expected from 
the fact that most loops accepted in the short-loop algo- 
rithm are small, while large loops that wind around the 
boundary of the lattice is needed to change the total po- 
larization. These are not suppressed in the directed-loop 
algorithm, thus leading to better performance. 

The Eqs. and @ form several coupled sets which in 
most cases are under-determined. There are thus many 
solutions for the out-edge selection probabilities P. Some 
general solutions and analysis of their efficiency for dif- 
ferent quantum systems were reported in Ref. Here 
we employ the solution B in Ref. to the eight-vertex 
model, but solutions for higher-vertex models are not 
hard to find as well. The allowed vertices for the eight- 
vertex model and their statistical weights are shown in 
Fig. [3 To shorten notation, we consider the so-called 
symmetric case: the statistical weights, a, b, c, and d, 
of the allowed states are assumed to be invariant under 
the simultaneous reversal of all arrows. The generaliza- 
tion of the algorithm to the non-symmetric case can be 
performed easily. 

Let Wi, . . . , Wi be the vertex weights a, b, c, d of the 
eight-vertex model ordered so that W\ > Wi > W 3 > 
W4,. Then the probability for picking the out-edge on a 
vertex with weight Wi resulting in a new vertex weight 
Wj after flipping arrows is Uj/Wi, where ty = tji and 
the non-zero entries of the 4x4 matrix t are 



tu = (Wi + W 2 - W 3 - W 4 )/2, 
tia = (Wi - W 2 + W 3 - Wa)/2, 
t 23 = (-W 1 + W 2 + W 3 + W 4 )/2, 



(3) 
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FIG. 4: The domain wall boundary conditions. Shown is an 
N x N lattice. The total number of vertices is iV 2 . The x 
and y coordinates take integer values at the midpoints of the 
horizontal edges. 



FIG. 3: Integrated autocorrelation times for number of 
c-type vertices (upper panel) and the total polarization in 
the y-direction (lower panel) for the directed-loop algorithm 
(open symbols) and the short-loop Barkema-Newman algo- 
rithm (closed symbols). The data shown is for the symmetric 
six-vertex model on an W x JV square lattice with PBC and 
vertex weights a = 6 = 2 and c = 1. 

tu = Wt, 

when W\ — Wi — W3 — W4 < 0. Otherwise one needs to 
include bounces in which the out-edge coincides with the 
in-edge. In this case a solution can be chosen as follows: 

t n = W1-W2-W3- W 4 , 

tij = t j i=W j , j = 2, 3, 4, (4) 

tij = 0, otherwise. 

The directed-loop algorithm satisfies ergodicity as any 
configuration can be obtained from another configuration 
by flipping spins along a finite number of (possibly over- 
lapping) loops. This process is exactly the directed-loop 
update, thus ergodicity follows. 

The algorithm presented here has many similar fea- 
tures to the loop algorithm ||. The loop algorithm 
breakup rules for the symmetric eight-vertex model can 
be chosen identical to Eq. as can be seen from 

Ref. 0, whenever the weights are such that no bounces 
are needed in the directed-loop algorithm. However in 
parameter regimes where bounces are needed, the related 
feature in the loop algorithm is to "freeze" independent 
loops together. Bounces and "freezing" of loops are very 
different in how they act to change the configuration. 
While bounces is a local resistance to changing a vertex, 
"freezing" causes big non-local changes of the vertex con- 
figuration. There are also other differences: For general 
vertex models the set of non-freezing/bouncefree solu- 
tions is always smaller for the loop algorithm than for 



the directed-loop algorithm. This allows the directed- 
loop algorithm to be efficient in a larger region of param- 
eter space than the loop algorithm. In particular this 
applies to the asymmetric eight-vertex model. 

Note that the need for bounces is generally not so cru- 
cial for higher-vertex models with many weights of the 
same magnitude, thus we expect that the directed-loop 
algorithm should work well in simulating these. Note 
also that an algorithm based on the directed-loop idea 
was recently demonstrated to be effective in simulating 
classical integer- valued link-current models |10| . 



III. SIX- VERTEX MODEL WITH THE DWBC 

The six-vertex model with the DWBC was introduced 
in Ref. in connection with the calculation of the cor- 
relation functions for exactly solvable 1 + 1 dimensional 
models . Here we recall the definition of the model in 
brief, referring for further details to the Ref. 01 where 
a more detailed description of the model and a compre- 
hensive list of references are given. 

The model is defined on an N x N square lattice; the 
thermodynamic limit corresponds to N — > 00. There are 
six possible states at each vertex: one should set d = 
in the eight-vertex model defined above, Fig. [21 The 
model is symmetric: the statistical weights, a, 6, and c, 
of the allowed states are assumed to be invariant under 
the simultaneous reversal of all arrows. Hence, the model 
is characterized by only two parameters, which can be 
taken to be a/c and b/c. We set c = 1 henceforth. 

The DWBC imply that all arrows on the top and bot- 
tom of the lattice are pointing inward, while all arrows 
on the left and right boundaries are pointing outward, 
Fig.H 

To investigate the spatially inhomogeneous behavior of 
this model we focus on the polarization, xn(x, y) |l2lll3| . 
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FIG. 5: The phase diagram of the six- vertex model in terms 
of the weights a and b. One has A > 1 in the regions I and 
II, — 1 < A < 1 in the region III, and A < — 1 in the region 
IV. The dotted quartercircle corresponds to A = 0. 



one finds by setting d = W4 = and c = 1 in Eqs. 10, 
that bounces are only necessary when a+b < 1 or |a— b\ > 
1. For the boundary vertices the loop is not allowed to 
exit on the boundary edges, because the arrows on these 
edges are fixed by the boundary conditions. This leads 
to more restricted equation sets (many W's are equal to 
zero) for the boundary vertices and generally requires the 
inclusion of bounce processes. 

Another important point should be mentioned is that 
the DWBC do not violate the ergodicity of the algorithm 
even though loops which wind around the boundaries 
are excluded. These winding loops are needed in order 
to change the net polarization in the x- or y-direction. 
However, one can verify that the boundary conditions 
restricts the net polarization in both these direction to 
be zero for any configuration, so winding loops are not 
necessary to sample the full configuration space allowed 
by the boundary conditions. 



which is the ensemble average of the arrow direction on 
the edge with coordinates (x, y) on the N x N lattice. 
The coordinate system used is shown in Fig. 0] Due to 
the symmetry of the model it is sufficient to consider the 
polarization of the horizontal arrows only. The value +1 
(— 1) is assigned to an arrow pointing to the right (left) 
and the ensemble average is assumed to be normalized 
by dividing by the partition function. Therefore, Xn lies 
between —1 and 1. 

Obviously, xn is independent of the coordinates of the 
edge in case of PBC. For these boundary conditions xn 
is known in the thermodynamic limit, and exhibits fer- 
roelectric order, antiferroelectric (AF) order or no order, 
depending on the position on the (a, b) plane. Thus, three 
phases exist in the six- vertex model with PBC: ferroelec- 
tric, antiferroelectric, and disordered phase. In Fig. [S] 
the phase diagram on the (a, b) plane for the model with 
PBC is plotted (cf., Fig. 8.5 of Ref. @). 

Introduce a parameter A by the formula 

lab 

The case A > 1 (regions I and II in Fig. [SJ corresponds 
to the ferroelectric phase; the case — 1 < A < 1 (region 
III in Fig. [SJ) to the disordered phase; the case A < — 1 
(region IV in Fig. [5j to the AF phase. 

Fig. [3] may be considered as the phase diagram for 
the model with the DWBC, in the sense that the free 
energy takes a different analytic form in the regions I 
through IV (see Ref. 0] for details). But, in case of the 
DWBC the polarization xn depends on the position of 
the edge. In the next section we show numerical results 
for the polarization xn(x,u) of the horizontal arrows as 
the parameters a and b are varied. 

Making use of the directed-loop algorithm described 
in the previous section for simulation of the model with 
the DWBC one should treat vertices belonging to the 
boundary and the bulk vertices differently. In the bulk 



IV. RESULTS 

In this section we present the results of the simulations 
for the polarization xn(%, y) in the disordered, antiferro- 
electric and ferroelectric phases. 

(i) Disordered phase: — 1 < A < 1. First consider the 
particular case A = (dotted quartercircle in Fig.0). An 
exact expression for XN(x,y) in this case was obtained 
by Kapitonov and Pronko [T]| recently. To check our al- 
gorithm we have compared results for the polarization at 
the point a = b = 1 /V2 with the exact results of Ref. . 
The comparison can be seen in Fig. El where the polar- 
ization is shown as a function of x for different values 
of y and system sizes, N. One can clearly see that the 
boundary values of the polarization, ±1, extends a finite 
distance into the bulk and this distance depends on y. 
The areas where the polarization stays at its boundary 
values are termed "frozen" regions. Going further into 
the bulk, there is a transition to a "disordered" region, 
where apart from small wiggles due to the finite system 
size, the polarization changes smoothly. It is interesting 
to note that there never is any extended regime where the 
polarization is zero, as is the case for PBC. The transi- 
tion between the "frozen" and "disordered" regions gets 
sharper as the system size is increased, as can be seen by 
comparing the two panels in Fig. 

It is convenient to visualize the behavior of the po- 
larization using greyscale plots, where greyvalues are as- 
signed to values of XN(x,y) and each point (x,y) cor- 
responds to a location of the midpoint of a horizon- 
tal edge following the layout described in Fig. For 
a = b = 1/V2 such a plot is shown in Fig. E{a). The 
four "frozen" corners are clearly apparent. In these re- 
gions, the vertices are all of the same type, and are, from 
upper left to bottom right, a\, bi, 62, respectively. 
To measure the area of the "frozen" regions, we define a 
threshold value e = 0.08, such that points (x, y) where 
I Xjv (a?, J/) I > 1 — e are judged to be in a "frozen" region. 




(x-l)/N 

FIG. 6: Polarization Xt*(x,y) as a function of x for different 
values of y. Vertex weights a — b = l/y/2. Results for two 
different system sizes are shown: N = 32 (upper panel) and 
N — 64 (lower panel). The filled symbols are Monte Carlo re- 
sults, while the crosses are exact results gotten from Ref. [l5|. 
The dotted lines are guides to the eye. 



Applying this we find that each "frozen" corner is 4.6% of 
the total area. This value changes relatively little chang- 
ing the value of e. 

Going away from the A = curve, let us follow along 
the diagonal, a = b, towards A = oo first, Fig. As 
the values of the vertex weights a and b increase, the 
area of the "frozen" regions decreases. We find that with 
e = 0.08 each frozen corner in (b) is 4.0% of the total 
area, and 2.8% in (c). For very large values of a = b, 
the polarization XN{x,y) increases linearly from —1 to 
1 as (x ~ 1)/N goes from to 1, independent of y, as 
can be seen in Fig. [7| (d) . This is consistent with what 
is expected from an ensemble of configurations with the 
smallest possible number of c-type vertices: Nl configu- 
rations each with a single c-type vertex on every row and 
column. 

Consider now a ^ b. Because of the symmetry of the 
phase diagram, Fig. |3J one can choose b > a without 
loss of generality. The weights of the vertices in the 
four "frozen" corners are no longer equal, and the "dis- 
ordered" region distorts into an oblong shape oriented 
along the diagonal with large corners of 62 and b\ vertices, 
see Fig. El The simulations for a = 1/4 and b = 
are shown in Fig. Eta)- The width of the oblong region 
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FIG. 7: Greyscale plot of the polarization xn(x, y) for N = 64 
in the disordered phase. Vertex weights are equal, a = b, and 
run through the values I/a/2, 1, 3, 100 for figures (a)-(d), 
respectively. The corresponding values of A are 0, 1/2, 17/18, 
1-5- 10" 5 . 
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FIG. 8: Greyscale plot of the polarization XN(x,y) for N = 
64. The weight a = 1 /4, while the weight b is chosen to be b = 
a/15/16 (A = 0, disordered phase) in figure (a) and b = 5/4 
(A = 1, the boundary between disordered and ferroelectric 
phases) in figure (b). 

shrinks as b increases keeping a fixed, a — 1/4, and be- 
comes very thin at the boundary to the ferroelectric re- 
gion, as can be seen in Fig. Efb). Along this boundary, 
b = a + 1, the width of the oblong region expands as a 
increases with N being constant. 

(ii) Antiferroelectric phase: A < — 1. The simulations 



FIG. 9: The two configurations having maximum number of 
the c-type vertices. These vertices are marked by filled circles. 
The size of the lattice is 4 x 4. 



in the AF phase are less efficient than in the disordered 
phase. This is partly due to the presence of the bounce 
processes also for bulk vertices, but another feature which 
makes the simulations difficult in this phase is the degen- 
eracy of the two types of AF orders. In the AF phase 
it becomes energetically favorable to have a maximum 
possible amount of c-type vertices, which is achieved by 
placing c-type vertices in a diamond placed in the center 
of the lattice. For an even N this diamond can be placed 
in two equivalent places differing only by one lattice spac- 
ing, as shown in Fig. [5] The Monte Carlo algorithm is 
however slow in tunnelling between these configurations, 
and this sets a limit to its performance. For odd N there 
is no such a degeneracy and the simulations are more effi- 
cient. Greyscale plots of the polarization for a = b = 1/2 
and a = b = 3/8 are shown in Fig. ^| We have plotted 
results for both even and odd N. 

One can see that the "disordered" region have a 
diamond-like shape, which is consistent with the dom- 
ination of the c-type vertices in the AF phase. As a = b 
decreases (A — » — oo), the shape of the "disordered" re- 
gion should converge to the one shown in Fig. El that is, 
the boundaries of the "disordered" region should become 
more and more straight. But, this convergence appears 
to be rather slow and it is not easy to see it from Fig. 1101 
What one can clearly see from Fig. is the difference 
between odd and even N. For odd N AF oscillations are 
clearly visible in the center of Figs. HOf cl and (d), while 
they are much weaker for even N, Figs. llOf a) and (b), 
reflecting the degeneracy mentioned above. These differ- 
ence between even and odd N can also be clearly seen 
from Fig. El For odd N AF oscillations are weaker at 
a = b = 1/2 than at a = b = 3/8. 

For a b greyscale plots are shown in Fig. m Here 
AF oscillations in the middle of the plot are visible for 
a = 1/4 and b = 1/2, Fig. Illf a'l. while they have almost 
vanished at the boundary between the AF and disordered 
phases, Fig. lriT b). 

(hi) Ferroelectric phase: A > 1. The behavior of the 
polarization in this phase is essentially the same as shown 
in Fig. |H1 Vertices of type b dominate completely in the 
region II of the phase plane Fig. |3J while in the region 
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FIG. 10: Greyscale plot of the polarization Xn(x,v) for two 
different system sizes: N = 32 in figures (a) and (b) and 
N = 33 in figures (c) and (d). Vertex weights are equal, 
a = b, and take the value 1/2 (A = — 1, the boundary between 
disordered and AF phases) for figures (a) and (c), and the 
value 3/8 (A = -23/9, AF phase) for figures (b) and (d). 



I of the phase plane the dominant vertices are those of 
type a. If one goes along the phase boundary, b = a + 1, 
towards a = oo, the widths of the "disordered" region is 
increased, as we have mentioned in the end of the part 
(i) of this Section. 

The exact expression is known 0] f° r the polariza- 
tion along the boundary xn{x, 1). Comparing our Monte 
Carlo data to this expression we find that in no cases is 
the absolute difference bigger than 0.016, which is com- 
parable to the statistical errors of our simulations. 



V. DISCUSSION 

We have considered the phase diagram of the model 
for the given N. Now, discuss the following problem: 
what happens with xn(x, y) in the thermodynamic limit, 
N — > oo? It is natural to expect that differences in the 
behavior of the polarization in the different parts of the 
phase plane, Fig [3] become more pronounced as N — ► oo. 
As one can see in Fig. the wiggles in the "disordered" 
region decrease with N increasing, and this is, indeed, the 
case for all the points (a, b) lying in the disordered phase 
(— 1 < A < 1, region III of the phase plane, Fig. and 
checked in our simulations. We expect that these wiggles, 
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FIG. 11: Greyscale plot of the polarization XN{x,y) for TV = 
32. The weight a = 1/4 while the weight b is chosen to be 
b = 1/2 (A = —11/4, AF phase) in figure (a) and 6 = 3/4 
(A = — 1, the boundary between disordered and AF phases) 
in figure (b). 




(x-l)fN 

FIG. 12: Boundary polarization xn(x,1) is shown for three 
system sizes, TV = 16, 32, and 64. Vertex weights a = 1/4 
and 6 = 3/4 (A = —1, the boundary between disordered and 
AF phases) . Note the steepening of the curve as TV increases. 



coming from the antiferroelectrically ordered configura- 
tions, should vanish completely in the thermodynamic 
limit for this phase. The next conjecture we want to 
make is on the behavior of the polarization along the 
boundary, xn{x, 1). It is known that for A = 0, as well 
as at the point a = b = 1, the boundary polarization be- 
comes the Heaviside step function in the thermodynamic 
limit [13, IT^. We conjecture that this is the case for the 
whole disordered phase; the position of the discontinuity 
will depend on the ratio between a and b. We present 
Fie. 1121 to support this conjecture. 

Furthermore, note that for a = b = l/v2 there is a 
mapping (see, e.g., Ref. |14j) of the six- vertex model with 
the DWBC onto the so-called model of domino tilings 
of the Aztec diamond. The thermodynamic behavior of 
the latter model was investigated in Refs. 17]. It shows 
the same features as in Fig. Ufa): the tilings are ordered 
(frozen) in the corners of the diamond, while going inside 
one falls into the "disordered" region. All these features 
were named the "Arctic Circle Theorem" , since the shape 
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FIG. 13: Polarization xn(x,d = TV/2 + 1) is shown for three 
system sizes, TV = 16, 32, and 64. Vertex weights a — 1/4 
and 6 = 5/4 (A = 1, the boundary between disordered and 
ferroelectric phases). Note the steepening of the curve as TV 
increases. 



of the boundary between the "frozen" and "disordered" 
regions is circular. The transition between "frozen" and 
"disordered" regions is step-like, with the height of the 
step function depending on the coordinates x and y. 

We expect the analogue of the Arctic Circle Theorem 
to take place for the whole disordered phase, — 1 < A < 
1: there should be the "frozen" regions, "disordered" 
region, and a sharp transition between them. We ex- 
pect also that the profile of the boundary between the 
"frozen" and "disordered" regions is circular for a = b, 
even though there is no obvious symmetry protecting 
this statement. Note that the very "smeared" profile 
in Fig. UKd) does not contradict our hypotheses because 
AT = 64 is relatively small compared to the values of the 
vertex weights a and 6, and is thus far from the thermo- 
dynamic limit for this point of the phase diagram. 

For the ferroelectric phase, A > 1, the greyscale plot 
Fig. ISfb) together with the scans shown in Fig. ED leads 
to the natural conjecture: in the whole region II of the 
phase plane, Fig.[5I a sharp discontinuity from a "frozen" 
domain with bi vertices to the one with &2-vertices takes 
place in the thermodynamical limit. In the region I the 
behavior is essentially the same, one should simply use 
a-type vertices instead of the 6-type. 

To this end, consider the antiferroelectric phase, A < 
— 1. We expect the step-like behavior of the boundary po- 
larization, \n{x,1) in this phase in the thermodynamic 
limit, as well as the existence of the "frozen" regions in 
the corners. Our statements on the behavior of the po- 
larization deep inside the lattice are more speculative. 
For a = b and even N the height of the AF oscillations 
decreases, while for odd N these oscillations seem to be 
non- vanishing in the large N limit, see Fig. [21 Our be- 
lief is that there is a finite region with AF order for odd 
N, as N — > oo, while for even N the polarization exhibits 
no such an order. 

Finally, we would like to stress that the directed-loop 
algorithm can also be applied to study the six-vertex 
model with any boundary conditions, and the higher- 
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FIG. 14: Polarization Xn{ x iV) along lines of constant y, 
where y = N/2 + 1, y = (N + l)/2) for even and odd N re- 
spectively, is shown for N = 8, 16, and 32 (upper panel) and 
N = 9, 17, and 33 (lower panel). Vertex weights a = 6 = 3/8 
(A = -29/9, AF phase). 



vertex models. These could help in solving the prob- 
lems for which the analytical methods are difficult to ap- 
ply. For example, the six-vertex model with any bound- 
ary conditions can be considered as a model for a de- 
scription of interface roughening of a crystal surface • 
An important point in studies Refs. |l8j is the existence 
of exact analytical results for the six-vertex model with 
PBC 0,0- Therefore, numerical data referring to other 
boundary conditions than PBC could give a new insight 
for these studies. 
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